Dynamic critical behavior of the Chayes— Machta— Swendsen— Wang algorithm 
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We study the dynamic critical behavior of the Chayes-Machta dynamics for the Fortuin-Kasteleyn 
random-cluster model, which generalizes the Swendsen- Wang dynamics for the q-state Potts model 
to noninteger q, in two and three spatial dimensions, by Monte Carlo simulation. We show that the 
Li-Sokal bound z > ot/v is close to but probably not sharp in d — 2, and is far from sharp in d = 3, 
for all q. The conjecture z > /3/u is false (for some values of q) in both d = 2 and d = 3. 
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Monte Carlo simulations in statistical mechanics |l| 
and quantum field theory typically suffer from criti- 
cal slowing- down 0,01 : the autocorrelation (relaxation) 
time r diverges as the critical point is approached, most 
often like r ~ £ z , where £ is the spatial correlation length 
and z is a dynamic critical exponent. For local algo- 
rithms, one usually has z ss 2. This effect severely limits 
the efficiency of Monte Carlo studies of critical phenom- 
ena in statistical mechanics and of the continuum limit 
in quantum field theory. 

An important advance was made in 1987 with the in- 
vention of the Swendsen-Wang (SW) cluster algorithm 
5 1 for simulating the g-state ferromagnetic Potts model 
6 ,0] at positive integer q. The SW algorithm is based on 
passing back and forth between the Potts spin represen- 
tation and the Fortuin-Kasteleyn (FK) bond represen- 
tation [1, 0| . This algorithm does not eliminate critical 
slowing-down, but it radically reduces it compared to lo- 
cal algorithms. Much effort has therefore been devoted, 
for both theoretical and practical reasons, to understand- 
ing the dynamic critical behavior of the SW algorithm as 
a function of the spatial dimension d and the number q of 
Potts spin states [10| . Unfortunately, it is very difficult to 
develop a physical understanding from the small number 
of "data points" at our disposal: second-order transitions 
occur only for (d,q) = (2,2), (2,3), (2,4), (3,2) and (4,2) 
El- 

A further advance was made in 1998 by Chayes and 
Machta (CM) [12], who devised a cluster algorithm for 
simulating the FK random-cluster model [H, [l3| — which 
provides a natural extension of the Potts model to non- 
integer q — at any real q > 1. The CM algorithm gener- 
alizes the SW algorithm and in fact reduces to (a slight 
variant of) it when q is an integer. By using the CM 
algorithm, we can study the dynamic critical behavior 
of the SW-CM dynamic universality class as a func- 
tion of the continuous variable q throughout the range 
1 < 9 < 9e(£), where q c (C) is the maximum q for which 
the transition is second-order on the lattice L i.14]. This 
vastly enhances our ability to make theoretical sense of 



the numerical results. 

In this Letter we report detailed measurements of 
the dynamic critical behavior of the CM algorithm 
for two-dimensional random-cluster models with 1 < 
q < 4 15] and for three-dimensional models with q = 
1.5, 1.8, 2, 2.2 0. Among other things, we find strong 
evidence against the conjecture z > (3/v recently pro- 
posed by two of us [n| , which had seemed plausible from 
the data for integer q. 

The FK random- cluster model with parameter q > is 
defined on any finite graph G — (V, E) by the partition 
function 
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where A is the set of "occupied bonds" and k(A) is 
the number of connected components ( "clusters" ) in the 
graph (V,A); here {v e } are nonnegative edge weights. 
For q = 1 this reduces to independent bond percolation 
[l|| with occupation probabilities p e = v e /(l + v e ); for 
integer q > 1 it provides a graphical representation of the 
g-state ferromagnetic Potts model with nearest-neighbor 
couplings { J e }, where v e = e^ J " — 1. 

It is convenient to consider a generalized random- 
cluster (RC) model 3, Tj^ 
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where Hi , . . . , are the connected components of the 
graph (V,A), and {W(H)} are nonnegative weights as- 
sociated to the connected subgraphs H of G. The model 
© reduces to the FK model (p} if W(H) = q for all 
H; other special cases include an FK representation for 
the Potts model in a magnetic field [2(| and various loop 
models [H. 

Now let m be a positive integer, and let us decom- 
pose each weight W(H) into m nonnegative pieces, any 
way we like: W(H) = YZ=i w <*{H). The first step of 
our generalized Chayes-Machta algorithm, given a bond 
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configuration A, is to choose, independently for each con- 
nected component Hi, a "color" a S {1, ...,m} with 
probabilities W Q (/fj)/W(ifi); this color is then assigned 
to all the vertices of Hi. The vertex set V is thus parti- 
tioned as V — UI!Li ^o* ^ ^ s n0 * nar d to see that, condi- 
tioning on this decomposition, the bond configuration is 
nothing other than a generalized RC model with weights 
{W a (H)} on the induced subgraph G[Vy, independently 
for each a. 

We now have the right to update these generalized RC 
models by any valid Monte Carlo algorithm. One valid 
update is "do nothing" ; this corres pon ds to the "inac- 
tive" colors of Chayes and Machta [12J. Of course, we 
must also include at least one nontrivial update. The 
basic idea is to have at least one color for which the 
weights W a (H) are "easy" to simulate. In particular, 
when W(H) = q for all H (the standard FK random- 
cluster model), we can take W a (H) = 1 for one or more 
colors a (the so-called "active" colors); the correspond- 
ing model on G[V a ] is then independent bond percola- 
tion, which can be trivially updated. Since we must have 
W a (H) < W(H), this works whenever q > 1. More gen- 
erally, if q > k, then we can have k active colors. If q is 
an integer and we take k = q, we recover the standard 
SW algorithm. 

Wc used the CM algorithm to simulate the random- 
cluster model in dimensions d = 2, 3 on hypercubic lat- 
tices of size L d with periodic boundary conditions. We 
measured the "energy-like" observables Af = # of occu- 
pied bonds and £' = # of nearest-neighbor pairs be- 
longing to the same cluster; the cluster-size moments 
Sm — Yl |C| m f° r various values of m, where |C| is the 
number of sites in the cluster C; and the size Ci of the 
ith-largest cluster for i = 1,2,3. 

For any observable O, let po(i) be its normalized au- 
tocorrelation function. Then define the exponential au- 
tocorrelation time 
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Figure 1: Integrated autocorrelation times Tj llti £/ versus lat- 
tice size L for the critical two-dimensional random-cluster 
model at 1.25 < q < 4, simulated using the Chayes-Machta 
algorithm with k — 1. 
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T cxp .o = limsup- 

t^±oo - log I po(t) I 
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Table I: Dynamic critical exponents Zi nt ,£' for two- 
dimensional random-cluster model as a function of q, with 
preferred fit and minimum L value used in the fit. Error bars 
are one standard deviation, statistical error only. The exact 
exponents a/u and pjv are shown for comparison [23j ]. 
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Typically all observables O (except those that, for sym- 
metry reasons, are "orthogonal" to the slowest mode) 
have the same value T CX p,e> = T cxp . However, they may 
have very different amplitudes of "overlap" with this 
slowest mode; in particular, they may have very differ- 
ent values of the integrated autocorrelation time, which 
controls the efficiency of Monte Carlo simulations 
We define dynamic critical exponents z oxp and z- mtt o by 
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On a finite lattice at 



We began by performing simulations on the square 
lattice (d = 2) at the exact critical point v c (q) = ^Jq 
2l| for 1.25 < q < 4 in steps of 0.25 and lattice sizes 
16 < L < 1024, using all positive integer values of 
k < q. We estimated the integrated autocorrelation times 
Tmt.o by the automatic windowing method described in 



criticality, £ can here be replaced by L. 



[10L [22j. The complete set of runs used approximately 
14.8 yr CPU time on a 1266 MHz Pentium III Tualatin 
processor. 

The autocorrelation functions of AT, £' and S m are 
in all cases very close to a pure exponential. In Fig. [T] 
we plot Tints' (for k — 1) versus L, and in Table [J we 
report the estimated dynamic critical exponents Zi n t,f- 
Our data also show that, as expected, the exponents are 
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Figure 2: Integrated autocorrelation times Ti nti £< divided by 
specific heat Ch, versus lattice size L, for the critical two- 
dimensional random-cluster model at 1.25 < q < 4, simulated 
using the Chayes-Machta algorithm with k = 1. 
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Figure 3: fcr inti g/ versus lattice size L for Chayes-Machta 
simulations of the three-dimensional random-cluster model 
with q — 1.5, 1.8, 2.2 at three near-critical temperatures, tak- 
ing k = Iq]- Data for q = 2, k = 2 at approximate criticality 
are from [10]. 



independent of fc, and we have roughly r oc 1/fc. 

Since the Li-Sokal bound Tin^Af, Ti n t,£' > const x Ch 
and hence z- lnt j^, z- ultt £' > a/v, originally proven [l7| for 
the Swendsen-Wang algorithm, can also be proven [l6| 
for the Chayes-Machta algorithm (at least for Af), it 
is of interest to analyze its possible sharpness [lH- In 
Fig. [2] we plot Tint,£'/C# versus L, in an attempt to de- 
termine whether this ratio is bounded or not as L — > oo. 
The results are far from clear, but our best guess is that 
T int,f ' /Ch diverges as L — > oo, either as a small power or 
as a logarithm. However, the precise behavior needs to 
be explored by simulations at larger L. 

On the other hand, Ossola and Sokal [l(| recently con- 
jectured, on the basis of the "data points" (d, q) — (2, 2), 
(2,3), (2,4), (3,2) and (4,2), that z ixA ,£> > (3/v; and 
they even speculated that we might have the equality 
Zint,£> — max(a/v, (l/v). The data for noninteger q now 
shed light on this conjecture: for q = 1.25, 1.5 there is 
modest evidence (and for q = 1.75 there is weak evi- 
dence) that z- mt £ i < (3/v, i.e. that even the weak form of 
the Ossola-Sokal conjecture is false. 

We next performed simulations on the simple-cubic lat- 
tice (d = 3) for q = 1.5, 1.8, 2.2 (see also [10] for q = 2) 
and lattice sizes 4 < L < 256, using k = the largest 
integer < q. We located the critical point by a finitc- 
size-scaling analysis using the ratio R — (S4)/ as in 
[25( 1 . The complete set of runs used approximately 21.5 yr 
CPU time on a 3.2 GHz Xeon EM64T processor. 

The autocorrelation functions of N, £' and S m are 
again very close to a pure exponential. In Fig. [3] we plot 
kT m t,£' versus L (multiplying by k makes the results for 
different q comparable) for three temperatures very near 
criticality. In Table [ill we report the estimated dynamic 
critical exponents Zmt.f and static critical exponents a/v 
and /3/v. In Fig.2]we plot fcri nt ,£' /Ch versus L. It seems 
clear that, for all four values of q, the Li-Sokal bound is 
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Table II: Dynamic critical exponents «i n t,£' an d static ex- 
ponents a/v and P/v for three-dimensional random-cluster 
model. For q = 2, dynamic data are from [l(i| ] and static 
exponents are from [29 ] . 



far from sharp. On the other hand, from Table ITll it seems 
clear that for q — 1.5, 1.8 we have the strict inequality 
Zint,£' < ft/v, once again ruling out the Ossola-Sokal 
conjecture even in its weak form. 

The dynamic critical behavior of the SW-CM dynamic 
universality class in dimension d = 3 therefore remains 
a mystery. Clearly, some new physical principle, beyond 
the slow equilibration of the energy embodied in the Li- 
Sokal bound [17[, needs to be discovered. 

One clue might be provided by our analysis [3] of the 
CM algorithm on the comple te g raph (mean-field limit), 
generalizing the analysis in [27| of the SW algorithm. 
Taking k = 1 and defining a "magnetization" m to be the 
fraction of sites in the largest cluster, we obtain for 1 < 
q < 2 the approximate difference equation (generalizing 
Beq. (10)]) 
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where ml is the value of m after a sweep in which the 
active group contains the largest cluster, and t is the 
deviation from the critical temperature. Clearly q = 2 is 
a special case because the coefficient of the linear term 
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Figure 4: fci"int,£' /Ch versus lattice size L for Chayes- 
Machta simulations of the near-critical three-dimensional 
random-cluster model at q = 1.5, 1.8, 2, 2.2, with k = [q\. 
Data for q — 2 at approximate criticality are from [Tot ] . 



equals 1: we have (3 = 1/2 and z = 1, and it is clear 
from the derivations 2, 27 1 that z is actually For 
1 < q < 2, by contrast, both the statics and dynamics are 
in the percolation universality class with (3—1 and z = 0: 
small perturbations from equilibrium relax exponentially 
with a finite autocorrelation time T cxPim = qj log[g/ {2q — 



2)] that diverges as q f 2. We conjecture that a similar 
behavior holds above the upper critical dimension, which 
for q < 2 is presumably <i = 6. Our numerical data 
[l6| confirm the behavior z = for 1 < q < 2 with 
t cx 1/(2 — g) as <? f 2, but not the predicted amplitude. 

Details of these simulations and their data analysis will 
be reported separately [H, [l6| ■ 
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